Two-dimensional nonlocal vortices, multipole solitons and azimuthons in dipolar 

Bose-Einstein condensates. 
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We have performed numerical analysis of the two-dimensional (2D) soliton solutions in Bose- 
Einstein condensates with nonlocal dipole-dipole interactions. For the modified 2D Gross-Pitaevski 
equation with nonlocal and attractive local terms, we have found numerically different types of 
nonlinear localized structures such as fundamental solitons, radially symmetric vortices, nonrotat- 
ing multisolitons (dipoles and quadrupoles), and rotating multisolitons (azimuthons). By direct 
numerical simulations we show that these structures can be made stable. 
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I. INTRODUCTION 

The recent first experimental realization of a degen- 
erate dipolar atom gas where a Bose-Einstein con- 
densate (BEC) of 52 Cr atoms has been observed, and 
optimistic perspectives in creating a degenerate gas of 
polar molecules have stimulated a growing interest in 
the study of BEC with nonlocal dipole-dipole interac- 
tions [1, H, @ . Dipole-dipole forces are anisotropic and 
long range, so that the inter-particle interaction becomes 
essentially nonlocal. 

Nonlocal nonlinearity naturally arises in many areas 
of nonlinear physics and plays a crucial role in the dy- 
namics of nonlinear coherent structures. In particular, 
a rigorous proof of absence of collapse in arbitrary spa- 
tial dimensions during the wave-packet propagation de- 
scribed by the nonlocal nonlinear Schrodinger equation 
(NLSE) with sufficiently general symmetric response ker- 
nel has been presented in Refs. Stable vortex 
dipole 0, HO, [UGl and azimuthon [HOSEi solitons 
in media with nonlocal nonlinear response were theoret- 
ically predicted. Finally, nonlocality induces attraction 
between solitons and allows for the formation of bound 
states of out-of-phase bright solitons [l5| and dark soli- 
tons [jj]. 

A very attractive feature of BEC with dipole-dipole 
interactions is that the interplay between the nonlocal 
interaction, which is only partially attractive and may 
be tuned by means of rotating orienting fields [l^ . and 
the usual local short-range contact forces, leads to the 
possibility of experimental realization of highly control- 
lable and stable solitary structures in BEC Q. 

Recently, Pedro and Santos [3| have studied the physics 
of bright solitons in two-dimensional (2D) dipolar Bose- 
Einstein condensates with repulsive short-range interac- 
tioins. Using the reduction procedure, they have ob- 
tained 2D modified Gross-Pitaevski equation with the 
nonlocal term describing dipole-dipole interaction and 
showed that the existence of stable 2D solitary waves 



is possible. 

In this paper, using 2D model suggested by Pedri and 
Santos [H, we study 2D solitary waves (SWs) in BEC 
with attractive short-range and nonlocal dipole-dipole in- 
teractions. As is known, the collapse of BECs at some 
critical number of atoms is the main consequence of the 
attractive nonlinearity The presence of nonlocal 

interaction, however, significantly changes the situation 
and leads to stable localized states. We present different 
types of SWs (fundamental solitons, vortices, nonrotat- 
ing and rotating multisolitons) and by direct numerical 
simulations show that these localized structures can be 
made stable. 



II. MODEL AND BASIC EQUATIONS 

A dipolar BEC, consisting of N particles with the 
dipole moment d oriented along the z-axis, at sufficiently 
low temperatures is described by a NLSE with nonlocal 
nonlinearity 
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where ty(r,t) is the condensate wave function normal- 
ized to the total number of particles: / ^(r)! 2 dr = Af. 
The coupling constant g corresponds to the local con- 
tact interaction and g = 4irh 2 a/m, where a is the s-wave 
scattering length. In the following, we consider a < 0, i.e. 
attractive short-range interactions. An external trapping 
potential is assumed to be of the form U(r) = muj 2 z 2 /2, 
with no trapping in the xy-pl&ne. All dipoles are as- 
sumed to be oriented along the trap axis. The nonlo- 
cal potential is due to the dipole-dipole interaction, and 
the kernel V(r) is given by V(r) = <?d(l — 3cos 2 0)/r 3 , 
where g<i = afigd 2 /4/k, 9 is the angle between the vector 
r and the dipole axis, [io is the magnetic permeability of 
the vacuum, and — l/2<a<lisa tunable parameter 

ana. 
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Assuming the anzatz 'l'(r) = ip(r±)<po(z), where the 
function (po(z), describing the condensate in the direc- 
tion of the tight confinement, is the ground state of the 
ID harmonic oscillator in the z-dircction and normaliz- 
ing the length, time, and wave function by l z /^/2, l/ui z , 
and (7V//2) 1 / 2 , respectively (where l z = {h/mojz) 1 / 2 ), 
authors of Ref. [|[ , following the standard reduction pro- 
cedure, obtained the following 2D equation 

i^t = -A x iP + gipl.\iP\ 2 + l3 J R(r-r')\ip(r')\ 2 dA, 
with the two free dimensionless parameters 
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(3) 



The Fourier transform of the kernel R(r) in Eq. @ is 
R(k) = 2 - 3 V^fce fc2 erfc(fc) , (4) 

where erfc(x) is the complementary error function, so 
that 



R(t) = 



1 



(27T) 



e* r i?(k) dk. 



(5) 



In what follows, since Eq. ||3J) admits an additional rescal- 
ing, the parameter g has been fixed at g — =Fl, where 
the — (+) sign corresponds to attractive (repulsive) short- 
range interaction. 

Equation ([2|) conserves the 2D norm, N = J \ip\ 2 dxdy, 
and energy 
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III. MODULATIONAL INSTABILITY 

An important feature of the dipole-dipole interaction is 
that, due to the anisotropy, it is only partially attractive. 
Correspondingly, the spectrum R(k) of the response func- 
tion R(r) is not sign definite (note, in this connection, 
that an analysis of 2D soliton dynamics in the frame- 
work of Eq. © with somewhat resembling Eq. ([Jl, but 
positive definite, kernel R(k) was performed in Ref. d3|). 
Equation @ has a solution in the form of plane wave 



^0 = l^o I exp(ik • r — iujt) 



(7) 



provided ui = fc 2 - g(l + (3 J i?(r)dr)|* | 2 - The stability 
properties of the plane wave essentially depend on sign 
dcfiniteness of the spectrum of nonlinear response func- 
tion 0, [n| . On the other hand, modulational instability 
(MI) (instability of the plane wave with amplification of 
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FIG. 1: The growth rates 7 of the modulational instability for 
(a) attractive (g = —1) and (b) repulsive (g — 1) short-range 
interactions. Values of (3 (the ratio between dipole-dipole and 
short-range interaction) are shown near the curves. 



both sidebands) is often considered as a precursor for 
the formation of bright solitons. Considering perturbed 
plane wave solutions in the form 



ip = (|^o I + Sip) exp(ik ■ r - iuxt), 



where 



Sip = ip + e l 
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(8) 



(9) 



and linearizing Eq. ^ around ^ m ^tp, one can obtain 
the growth rate 7 of MI of homogeneous field (ko = 0) 
for the model Eq. © 



7 2 = -2\y \ 2 k 2 g[l + rJR(k)] 



(10) 



Instability occurs if j 2 > 0. In Figure [T] we show the 
dependence of the growth rate of MI on k for the cases 
of attractive (g = — 1) and repulsive (g = 1) short-range 
interactions. In the attractive case, the growth rate is 
equal to zero for < k < k cr , where k cr is some critical 
value depending on (3 (the ratio between dipole-dipole 
and short-range interaction), if (3 < 0.5 (i. e., in par- 
ticular, for all negative (3), so that long wave modes are 
stable. Optimal, i.e. corresponding to maximum of the 
growth rate, wave number k op t decreases with increasing 
(3. In the repulsive case, the growth rate of MI is equal 
to zero for all positive (3 and for (3 > —0.4. This is in 
agreement with results of Ref. [H , where bright solitons 
(for the repulsive case g = 1) were predicted only for 
negative (3 and \(3\ > 0.12. 



IV. NUMERICAL RESULTS 

We look for stationary solutions of Eq. with 
g = — 1 (attractive short-range interaction) in the form 
ip(x,y,t) = ^{x, y) exp(— i/xt), where /i is the chemical 
potential, so that ^ obeys the equation 



= -tf |V| 2 + p3 J R(r - r')|*(r')| 2 dr'\ , 

(11) 
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FIG. 2: Numerically found stationary localized nonrotating (a-c) and rotating (d-f) solutions of Eq. ([2]) with /3 = 2: (a) 
monopole with /i = —2\ (b) dipole with [i — —2; (c) quadrupole with fi = —2; (d) vortex with fi = —5; (e) azimuthon with 
fi = —5, p = 0.6 and two intensity peaks; (f) azimuthon with n — —5, p = 0.9 and four intensity peaks. 



where R(r) is determined by Eqs. (j4]) and ([5|). To solve 
numerically Eq. (jTTJ) , we impose periodic boundary con- 
ditions on Cartesian grid and use the relaxation tech- 
nique similar to one described in Ref. [H]]. We have not 
found any localized solutions with jj, > 0. Fundamental 
soliton solutions of Eq. (TlTj) with (5 < or (3 > f3 cr , where 
(3 cr ~ 2.1, turn out to be unstable with dN/dfi > 0. 
Thus, in what follows, we consider the region < (3 < (3 cr 
and specifically set f3 = 2. Choosing an appropriate ini- 
tial guess, one can find numerically with high accuracy 
(the norms of the residuals were less than 10~ 9 ) three dif- 
ferent classes of spatially localized solutions of Eqs. (fTTj) 
- the nonrotating (multi)solitons, the radially symmetric 
vortices, and the rotating multisolitons (azimuthons). 

The real (or containing only a constant complex fac- 
tor) function ^(x^y) corresponds to nonrotating solitary 
structures. Examples of such nonrotating (multi)solitons 
for Eq. pip , namely, a monopole, a dipole, and a 
quadrupole are presented in Figs.[2ja)-[2tc). The nonro- 
tating multipoles consist of several fundamental solitons 
(monopoles) with opposite phases. 

The second class of solutions, vortex solutions, are 
the solutions with the radially symmetric amplitude 
|\&(x, that vanishes at the center, and a rotating spi- 
ral phase in the form of a linear function of the polar an- 
gle 9, i.e. arg ^ = m6, where m is an integer. The index 
m (topological charge) stands for a phase twist around 
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FIG. 3: (a) Chemical potential /j, and (b) energy E versus nor- 
malized number of atoms N. Solid line: nonrotating dipoles 
(p = 0). Dotted line: vortices (p = 1). 



the intensity ring. The important integral of motion as- 
sociated with this type of solitary wave is the angular mo- 
mentum, which can be expressed through the vortex am- 
plitude and phase. The numerically found single-charged 
(to = 1) vortex solution of Eq. (fTTj) is shown in Fig.[^d). 

The third class of solutions, rotating multisolitons with 
the s pat ially modulated phase, were first introduced in 
Ref. |20| for models with local nonlinearity, where they 
were called azimuthons. The azimuthons can be viewed 
as an intermediate kind of solutions between the rotating 
radially symmetric vortices and nonrotating multisoli- 
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tons. Using variational analysis to describe azimuthons, 
the authors of Ref. [12] considered the following trial 
function in polar coordinates (r,0) 

*(r, 9) = r |m| $(r) (cos m6 + ip sinm0), (12) 

where $ is the real function, which vanish fast enough 
at infinity, m is an integer, and < p < 1. The case 
p = corresponds to the nonrotating multisolitons (e. g. 
to = 1 to a dipole, m — 2 to a quadrupole etc.), while the 
opposite case p = 1 corresponds to the radially symmetric 
vortices. The intermediate case < p < 1 corresponds to 
the azimuthons. In our case, the numerically found com- 
plex function ^(a;,y) with a spatially modulated phase 
corresponds to the azimuthons. We introduced the pa- 
rameter p (modulational depth), which is similar to the 
one in Eq. ijl2|L in the following way 

p = max|Im*|/max|Re\f f |. (13) 

For fixed chemical potential //, there is a family of az- 
imuthons with different p. Like the radially symmetric 
vortices, the azimuthons carry out the nonzero angular 
momentum. In Figures [21(e) and &(f) we demonstrate 
two numerically found examples of the azimuthons for 
the nonlocal model described by Eqs. Figure [3] 

shows the dependences of the chemical potential /i and 
energy E on the normalized number of atoms TV, for the 
dipoles (p = 0) and vortices (p = 1). 

We next addressed the stability of these localized so- 
lutions and study the evolution of the solitons in the 
presence of small initial perturbations. We have un- 
dertaken extensive numerical modeling of Eq. ([2]) ini- 
tialized with our computed solutions with added gaus- 
sian noise. The initial condition was taken in the form 
$(x,y)[l + e$(x, y)], where &(x, y) is the numerically 
calculated exact solution, $(x,y) is the white gaussian 
noise with variance a 2 = 1 and the parameter of per- 
turbation e = 0.005 -j- 0.1. In addition, azimuthal per- 
turbation of the form is sin 9 was taken for the vortices 
and azimuthons. Spatial discretization was based on the 
pseudospectral method. Under this, since the Fourier 
transform of the kernel R is known, the nonlocal term 
can be easily computed with the aid of the convolution 
theorem. Temporal i-discretization included the split- 
step scheme. As was said above, we consider the region 
< /3 < I3 cr . 

The fundamental solitons are stable for all /i < 0. 
These solitons have dN/dfi < so that the Vakhitov- 
Kolokolov stability criterion is met. 

Depending on the parameter [i, we observed three 
different scenarios of the nonrotating dipole evolution, 
which are presented in Fig. 2] (for e — 0.01). The first 
regime corresponds to the region ji cr < ji < 0, and for 
/3 = 2 we found ji cr ~ —0.2, which corresponds to the 
normalized number of atoms N cr ~ 15.5. If fi cr < ft, 
the initial dipole splits in two monopoles which move in 
the opposite directions without changing their shape and 
without radiation, i. e. the monopoles just go away at 
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FIG. 4: Left column: Splitting of the dipole with fj, = —0.2 
into two monopoles; middle column: stable dynamics of the 
dipole with /i = —3; right column: unstable dipole with \i = 
-4. 



infinity. This type of the evolution is shown in the left col- 
umn of Fig. [U Under this, the value SN = Ndi P — 2N mon , 
where N^ip and N mon are 2D norms for the dipole and 
monopole respectively, tends to almost zero as fi ap- 
proaches ju, cr . 

The second regime of the dipole evolution corresponds 
to the region fi t h < M < Mcr, where (for (3 = 2) [i t h ~ 
—3.1. The numerical simulations clearly show that in 
this range of the parameter fj, the dipoles are stable with 
respect to initial noisy perturbations and survive over 
huge times. In terms of the 2D norm (normalized number 
of atoms), the stability region is written as N cr < N < 
N t h, where Nth ~ 71. The stable dynamics of the dipole 
is illustrated in the middle column of Fig. [4] (for (3 = 2, 
[i = — 3 and e = 0.01). 

The further (after fith ~ —3.1) decreasing of the chemi- 
cal potential fj, (or, equivalently, increasing of the normal- 
ized number of atoms N) sharply shortens the times at 
which the dipole survives, and, the dipoles with \i < fi t h 
are unstable. The typical decay of the unstable dipole 
below the threshold value fi t h of the chemical potential 
is shown in the right column of Fig. [4] Thus, the stable 
dipoles exist only within a finite, rather narrow range of 
the normalized number of atoms N. 

A somewhat different behavior we observed for the vor- 
tices. The numerical simulations clearly show that the 
vortices with \x < fJ, cr , where \i cr is some critical value 
and for f3 — 2 we found fj, cr ^ — 1.4 (with corresponding 
N cr ~ 45), are stable with respect to small initial noisy 
and azimuthal perturbations up to the maximum times 
used (of the order of t = 1000). The vortices with \x > \i cr 
(i.e. N < N cr ) splits in two fundamental solitons moving 
in the opposite directions. These two different scenarios 
of the vortex evolution are illustrated in Fig. [5] Thus, the 
vortices can be made stable if the 2D norm (normalized 
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FIG. 5: (a) Splitting of the vortex with p = — 1; (b) stable 
dynamics of the vortex with p = —5. 

number of atoms) exceeds some critical value iV cr . 

We have not performed numerical analysis of the az- 
imuthon evolution for different A and arbitrary p because 
of the difficulties in finding azimuthon solutions with ar- 
bitrary p. Nevertheless, we can conclude that azimuthons 
with two intensity peaks and not too small p can be made 
stable if the 2D norm (normalized number of atoms) N 
exceeds some critical value depending on p. Splitting of 
the azimuthon with two intensity peaks and fi = — 1, 
p = 0.4, and stable dynamics of the azimuthon with 
fx = —5 and p = 0.6 are shown in Fig. [5] Numerically 
estimated rotational velocity of the stable azimuthon in 
Fig. Ob) is lu = 0.18 so that it survives over many dozens 
of rotational periods. 

Note, that one point should be emphasized. Strictly 
speaking, our direct numerical modelling can not give 
a rigorous proof of stability/instability of the multisoli- 
tons. First, in the the direct numerical experiments one 
can consider the evolution over finite times only. Second, 
the results are limited to the perturbation profile. A 
rigorous proof could, for instance, include a linear stabil- 
ity analysis with the corresponding eigenvalue problem. 
Nevertheless, from our numerical simulations of the dy- 
namics over finite, but large, times we can conclude that 
(in stable cases) the potential growth rates of unstable 
modes are very small. The structures (if stable) survive 
over huge times and hundreds of rotational periods, and 
from the practical point of view they can be regarded as 
stable. 



dipole-dipole and attractive short-range contact interac- 
tions and studied their stability. We have found numer- 
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FIG. 6: (a) Splitting of the azimuthon with two intensity 
peaks and \i = — 1, p = 0.4; (b) stable dynamics of the az- 
imuthon with two intensity peaks and p = —5, p = 0.6. 



ically three kinds of soliton families: nonrotating multi- 
pole solitons (fundamental one-hump soliton, dipole and 
quadrupole), radially symmetric vortices, and rotating 
multihump (with two and four intensity peaks) solitons 
with the spatially modulated phase (azimuthons). We 
have shown that stable solitons may exist only within 
a finite range of the ratio between dipole-dipole and 
short-range interactions (both of which are tunable). 
The anisotropy of the dipole-dipole interaction is crucial, 
since this leads to partially attractive nature of the inter- 
action. Sufficiently large dipolar interactions destabilize 
the SWs. By direct numerical simulations, we have found 
that dipole nonrotating solitons, vortices and two inten- 
sity peak azimuthons can be stable for some values of the 
chemical potential (or, equivalently, normalized number 
of atoms). 
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